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Abstract —We introduce a Prony-like method to recover a 
continuous domain 2-D piecewise smooth image from few of its 
Fourier samples. Assuming the discontinuity set of the image 
is localized to the zero level-set of a trigonometric polynomial, 
we show the Fourier transform coefficients of partial derivatives 
of the signal satisfy an annihilation relation. We present neces¬ 
sary and sufficient conditions for unique recovery of piecewise 
constant images using the above annihilation relation. We pose 
the recovery of the Fourier coefficients of the signal from the 
measurements as a convex matrix completion algorithm, which 
relies on the lifting of the Fourier data to a structured low- 
rank matrix; this approach jointly estimates the signal and the 
annihilating filter. Finally, we demonstrate our algorithm on 
the recovery of MRI phantoms from few low-resolution Fourier 
samples. 


I. Introduction 

The recovery of continuous domain parametric representa¬ 
tions from few measurements using harmonic retrieval/linear 
prediction has received considerable attention in signal pro¬ 
cessing since Prony’s seminal work 0 0 Extensive research 
has been devoted to the recovery of finite linear combination of 
exponentials with unknown continuous frequencies as well as 
linear combination of Diracs and non-uniform 1-D splines with 
unknown locations/knots 0 -j5|. Recently, convex algorithms 
that minimize atomic norm were also introduced to recover 
such continuous signals Q; these methods are off-the- 
grid continuous generalizations of compressed sensing theory, 
which can avoid discretization errors and hence are potentially 
more powerful. However, the direct extension of the above 
Prony-like and convex off-the-grid methods to 2-D piecewise 
smooth images is not straightforward. Specifically, the partial 
derivatives of piecewise smooth images can be thought of as 
a linear combination of a continuum of Diracs supported on 
the curves separating the regions, which current methods are 
not designed to handle. 

Recently, Pan et al. introduced a complex analytic 
signal model for continuous domain 2-D images, such that 
the complex derivatives of the image are supported on a curve. 
Under the assumption that the curve is the zero-set of a band- 
limited function, the authors show the Eourier transform of 
the complex derivative of such a signal is annihilated by 
convolution with the Eourier coefficients of the band-limited 
function. This property is used to extend the finite-rate-of- 


innovation model to this class of 2-D signals. This work 
has several limitations. One problem is that a complex analytic 
signal model is not realistic for natural images (e.g., the 
only real-valued analytic functions are constant functions). In 
addition, if we are only measuring a finite number of Eourier 
samples, one can choose analytic functions such that the all of 
these coefficients vanish; i.e. the recovery of the signal from 
few Eourier coefficients of an arbitrary signal in this class is 
ill-posed. 

In this work, we address the above limitations by proposing 
an alternative signal model based on the a more realistic 
class of piecewise smooth functions. We show that we can 
generalize the annihilation property in 0 to this class of 
functions, such that it is possible to recover an exact contin¬ 
uous domain representation of the edge set from few Eourier 
samples. We also determine necessary and sufficient conditions 
on the number of Eouier samples required for perfect recovery 
of the edge set in the caseof piecewise constant images. 

To recover the full signal, we propose a single-step convex 
algorithm which can be thought of as jointly estimating the 
edge set and image amplitudes. This is fundamentally different 
than the two stage approach proposed in which we 
also investigated for super-resolution MRI in ||^. Motivated 
by recently proposed algorithms for calibration-free parallel 
MRI recovery |T0| , G), our new approach is based on the 
observation that the ideal Eourier samples of the signal lift to a 
structured low-rank matrix, which allows us to pose the recov¬ 
ery as a low-rank structured matrix completion problem. We 
demonstrate our algorithm on the recovery of MR phantoms 
from low-resolution Eourier samples. 


II. Signal Model 
A. Piecewise smooth signals 

In this paper, we consider the general class of 2-D piecewise 
smooth functions: 


/(r) = Xni(r), Vr = (a;, y) € [0, l]^ (1) 




where xo denotes the characteristic function of the set U: 


Xn(r) = 


1 if r = (x, G U 
0 else. 


( 2 ) 


Here we assume each C [0,1]^ is a simply connected 
region with piecewise smooth boundary dVti. The functions gi 
in Q are smooth functions that vanish under of a collection of 
constant coefficient differential operators D = {Di^ 
within the region 

Dj gi(r) = 0,\lr j = l,..,N. ( 3 ) 

We now show that the above class of functions is fairly 
general and includes many well-understood image models by 
appropriately choosing the set of differential operators D. 

1) Piecewise constant images: We set D to 

D = V = 

Note that T>g = 0 if and only if g = Ci for q G C. Hence, Q 
reduces to the well-known piecewise constant image model: 

n 

/W = Vr = (a;,2/) G [0, l]^ (4) 

i=l 

This case will be the primary focus of this work due to its 
simplicity and provable guarantees. 

2) Piecewise analytic images: Choosing = dz = dx ^ 
jdy, then = 0 if and only if g is complex analytic. 
Hence this model is equivalent to the one proposed in Q. 
As described above, this signal model is not very realistic for 
natural images. 

3) Piecewise harmonic : Both the above cases consider only 
first-order differenial operators. One choice of a second-order 
differential operator is the Laplacian D = A = 5^^ -f dyy. 
Then = 0 if and only if g is harmonic. 

4) Piecewise linear images: If we consider all second or¬ 
der partial derivatives D = 9^^}, then = 0 if 

and only if g is linear, i.e. ^(r) = (a, r) b, for a G C^, 
6 G C, and so / has the expression 

n 

/(r) = bi) XQi (r), Vr = {x, y) G [0,1]^. (5) 

i=l 

5) Piecewise polynomial: Generalizing the above case we 
may consider all nth order partial derivatives D = {9^}|a|=n 
where a is a multi-index, then = 0 if and only if ^ is a 
polynomial of degree at most (n — 1). 

We will show that under certain assumptions on the edge set 
C = the Fourier transform of derivatives of a piece- 

wise smooth signal specified by 0 satisfies an annihilation 
property. This will enable us to recover an exact continuous 
domain representation the edge set C of a piecewise smooth 
signal from finitely many of its Fourier samples by solving a 
linear system. 

B. Trigonometric polynomials and curves 

Following ||^, we will assume the edge set C to be the 
zero-set of a band-limited periodic trigonometric polynomial 

^(r) = ^ c[k] (6) 

keA 

where c[k] G C and A is any finite subset of Z2; we call 
any function g described by 0 a trigonometric polynomial. 


and the zero-set C : {/i = 0} a trigonometric curve. We 
also define the degree of a trigonometric polynomial g to 
be the dimensions of the smallest rectangle that contains the 
frequency support set A, denoted as deg{g) = {K,L). For 
trigonometric polynomials g and u, we say u divides g or 
u \ gif g = where 7 is another trigonometric polynomial. 

Using elementary results from algebraic geometry, we may 
show there is a unique minimal degree trigonometric polyno¬ 
mial associated with any trigonometric curve C, which we call 
the minimal polynomial for C\ 

Proposition 1. For every trigonometric curve C there is a 
unique (up to scaling) trigonometric polynomial go with C : 
{go = 0} such that for any other trigonometric polynomial g 
with C \ {g = f)} we have deg (go) < deg{g) and go \ g. 

The following property of minimal polynomials is also 
important for our uniqueness results: 

Proposition 2. Let C be the zero set of a trigonometric 
polynomial with minimal polynomial go. Suppose v is a 
trigonometric polynomial such that u = ft and Vu = ft on 
C, then go \ u. In particular, Vgo = ft for at most finitely 
many points on C. 

HI. Annihilation property 

We now show that the Fourier transform of the partial 
derivatives of piecewise smooth signals 0 satisfy an anni¬ 
hilation property. 

1) First order partial derivative operators: First we con¬ 
sider the case of a single characteristic function xq. Note that 
since xn is non-smooth at the boundary, its derivatives are 
only defined in a distributional sense. Letting p denote any 
test function we have: 

{dxXQ,¥>) = -{xn,dx^) = - dx<fdr = - f '^dy ( 7 ) 

Ja Jan 

where the last step follows by Green’s theorem. Likewise, 

{dvXn,‘P)=f ‘Pdx (8) 

Jan 

Hence dxXn dyXn can be interpreted as a continuous 
stream of weighted Diracs supported on dLl. In particular, if 
is any smooth function that vanishes on dQ then 

V’ • dxXn = V’ • dyXn = 0 (9) 

where equality holds in the distributional sense. Assuming 
= g is Si trigonometric polynomial, taking Fourier trans¬ 
forms of (|^ yields the following annihilation relation: 

Proposition 3. Let f = Xn boundary dfl given by 
the trigonometric curve C \ {g = ft}. Let D be any first 
order differential operator. Then the Fourier transform of Df 
is annihilated by convolution with the Fourier coefficients 
c[k], k G A that is 

c[k] Df{(jj — 27rk) = 0, for all u; G (10) 

keA 








Due to the above property, we call /i an annihilating 
polynomial for Df. 

It is straightforward to extend the above proposition to 
piecewise constant functions / = provided /i = 0 

on the union of the boundaries C = Likewise, if 

f = g ' Qq where Dg = 0, then by the product rule 

Df = Dg ■ XQ+ 9 ■ Dxq. = 9 ■ Dxa 

which has support on dfl and so, fi ■ Df = 0, which implies 
holds for /, and similarly for the linear combination / = 

Yh=i 9i ■ Xcii, where Dgi = 0 for all i = 1 , 

2) Second order partial derivative operators: Now con¬ 
sider the case where D is any second-order differential opera¬ 
tor. Let f = P'Xn where Dg = 0. We now show that /i^ is an 
annihilating polynomial for Df, where g is any trigonometric 
polynomial that annihilates the partial derivatives of xn- 

Let d‘^ = 92^1 where di G {dx,dy},i = 1,2. By the 
product rule we have: 

d'^f = d‘^9 ■ Xn + dig ■ d2Xu + ^25 ' ^iXn + 9 ' d'^Xn- 
Since diXn and ^ 2 X 0 are annihilated by g, we have 

Again by the product rule 

g ■ d'^xn = d2{g ■ diXQ.) - 2/i • ^ 2 /^ • dixn = 0, 
which implies 

g-dg = xn-g-d^9 

and so by linearity 

g ■ Df = xn ■ g ■ Dg = 0. 

The above shows that g‘^ is always sufficient to annihilate 
Df, where D is second-order. However, using Prop. we 
may show that when g does not vanish on dQ, then g‘^ is 
also necessary for annhilation of Df, in the sense that if u is 
any other trig polynomial satisfying a • Df = 0, then g^ \ u, 
where g^ is the minimal polynomial for 9(2. If deg (go) = 
{K^L), this implies any annihilating polynomial a for Df 
has degiy) > {2K — 1, 21/ — 1). 

3) Partial derivative operators of arbitrary order: A simi¬ 
lar argument shows that when D is any nth order differential 
operator, and f = g ' Xn where Dg = 0, then 

g^'Df = 0. ( 11 ) 

This yields the following annihilation relation for higher-order 
differential operators: 

Proposition 4. Let D be any nth order differential operator. 
Let f = g ' Xo. ^ith Dg = 0 and 90 C {/i = 0} for some 
trigonometric polynomial g. Then the Fourier transform of 
Df is annihilated by convolution with the Fourier coefficients 
d[k], k G r o/ g^, that is 

d[k\ Dfiyj — 27rk) = 0, for all u; G (12) 

ker 

Likewise, by linearity, ( p^ is valid for linear combinations 
f = Y.igi' where Dgi = 0 and /i = 0 on [Jidfli. 


IV. Recovery from finite Fourier samples 

We now investigate necessary and sufficient conditions for 
the recovery of the filter coefficients describing the edge set 
from finitely many Fourier samples of the original signal /. 
For these results we restrict our attention to piecewise constant 
signals. 


A. Necessary conditions 

For a piecewise constant signal /, from the annihilation 
condition ( p^ we may form the linear system of equations: 

f SkeA ^ VIgF (13) 

\EkeA (27r[l - k]) = 0, 

where fx and fy may be computed from samples of / by 
Ix{g = -j^x ■ and Jy(u}) = -jujy ■ 7(w). Supposing 

the sampling grid (2 is a rectangular of dimensions {K',L'), 
and the minimal polynomial for C has degree (K^L), with 
coefficients c[k] supported in A, then we may form at most 
M = 2 • {K' — iX -f 1) • (!/' — L 1) valid equations from 
( p^ . Therefore to solve for the at most K • L unknowns c[k], 
k G A, we require at least M = K • L equations. This gives 
the following necessary condition for recovery of C\ 

Proposition 5. Let f be piecewise constant such that the edge 
set C has minimal polynomial g of degree (iX, L). A necessary 
condition to recover the edge set C from is to collect 
samples of f on a (iX',I/') rectangular grid such that 

2 • (iX' - iX + 1) . (L' - L + 1) > iX • L. 


To illustrate this bound, suppose the minimal polynomial 
has degree (iX, iX), and we take Fourier samples from a square 
region. Then this requires at least l.TliX x l.TliX Fourier 
samples to recover the edge set C. Our numerical experiments 
on simulated data (see Fig. [T]) indicate the above necessary 
condition might also be sufficient for unique recovery; that is, 
we hypothesize the minimal filter coefficients c[k] are the only 
non-trivial solution to the system of equations 

B. Sufficient conditions 

We now focus on sufficient conditions for the recovery of 
the edge set. Here we will use mA to denote a dilation of the 
set A by a factor of m: if A = : \k\ < iX, |/| < L}, 

then mA = {{k, 1) : \k\ <mK, |/| < mL}. 


Theorem 6. Let f = xn be the characteristic function of a 
simply connected region (2 with boundary 9(2 having minimal 
polynomial g with coefficients c[k],k G A. Then the c[k] can 
be uniquely recovered (up to scaling) as the only non-trivial 
solution to the equations 


EkeAc[k]L(27r[l-k]) = 0, 
EkeAc[k]7/(27r[l-k]) = 0, 


(14) 


The proof entails showing any other trigonometric poly¬ 
nomial p{r) having coefficients d[k\, k G A satisfying 
must vanish on 9(2, from which it then follows that 77 is 
a scalar multiple of the minimal polynomial g by degree 





(a) Original signal 



(b) Recovered fi 



(c) Edge set {ii = 0} 



(d) Singular values (log scale) 


Fig. 1. Exact recovery of edge set of a piecewise constant signal 
from the minimum necessary number of Fourier samples. The original 
piecewise constant signal is shown in (a), and was generated to have 
edge set C : {/i = 0} where the degree of ji is known to be (9,9). 
Here the minimal polynomial and edge set C is shown to be recovered 
from the system {T3} having dimensions 81 x 81 corresponding to 
the 15x15 necessary minimum number of Fourier samples predicted 
by Prop. 1^ The singular values of this system are plotted in (d), 
indicating the recovered is the only non-trivial solution. 


considerations. We also are able to show similar result holds 
for piecewise constant signals, provided the characteristic 
functions do not intersect: 


Theorem 7. Let /(r) = (r) be piecewise con¬ 

stant, where the boundaries dVti are described by non¬ 
intersecting trigonometric curves {p^i = 0}, where pi is the 
minimal polynomial for dVti. Then, the coefficients d[k],k G A, 
of p = pi' • • pn, and equivalently the edge set C = 
can be uniquely recovered (up to scaling) as the only non¬ 
trivial solution of 


EkeA<^M/r{27r[l-k]) = 0, 

EkeA«,(27r[l-k])=0, 


(15) 


Note that to form the ^uations in ( fT^ and requires 
access to Fourier samples /[k] for all k G 3A, which is greater 
than necessary number of samples given in Theorem We 
conjecture that the uniqueness results in Theorems [^and]^ can 
in fact be sharpened to the necessary number of samples, and 
extended to piecewise constant signals where the boundaries 
of the regions intersect. 


V. Recovery Algorithms 

Up to now we have only considered the problem of re¬ 
covering the edge set of a piecewise constant signal / = 
from finite Fourier samples. In analogy with 


Prony’s method, once the edge set is determined it is theoreti¬ 
cally possible to recover the signal amplitudes by substuting 
/ back into and solving a full rank system. However, 
this is not feasible in practice since it requires factoring a 
high degree multivariate polynomial into its irreducible factors. 
Instead we pursue approaches that allow us to pose the 
recovery as the solution to a convex optimization problem. 


A. Curve-aware recovery 

Supposing we have access to an annihilating polynomial 
p of the signal (equivalently, the edge set C : {/i = 0}), one 
approach is to pose the recovery as the weighted total variation 
minimization problem: 


/ = argmin / \p(r) V^(r)| dr subject to ^[k] = /[k], Vk G F 

9 J 

( 16 ) 

Here, since p = 0 on the edge set C, the gradient of the image 
not penalized along the curve, which allows for the recovery 
of an image with sharp edges along C. 

A version of this approach was investigated in an earlier 
work for the super-resolution recovery of MR signals from 
few Fourier samples j^. However, we found this scheme to 
have certain drawbacks, namely that it requires discretization 
onto a spacial grid, the optimization of many parameters, and 
is very sensitive to the estimate p. Hence we consider an 
alternate approach which does not rely on an explicit estimate 
of the annihilating polynomial p, but instead jointly recovers 
the image and the annihilating polynomial in a single stage 
algorithm. 


B. Low-rank recovery 

The annhilation equations specified by ( pA] ) can be repre¬ 
sented in the matrix form as 


T.[/] 

T,[/] 


d = 0 


T[/] 


(17) 


where d is a vectorized version of the Fourier coefficients 
(i[k],k G A, and and are block Toeplitz matrices 
corresponding to the 2-D convolution of d[k] with the discrete 
samples of jujx f{^) and jujy /(cj), respectively, for lv = 27rl, 
1 G U. Specifically, if the coefficient support set A has 
dimensions KxL, each row of (/) is the vectorized version 
of an K X L patch of juj^ f{^), and likewise for Ty{f). The 
number of rows is equal to the number of distinct patches, 
which correspond to the number of annihilation equations. 

Note that for a given piecewise constant signal, a priori we 
do not know the degree of the minimal polynomial describing 
the edge set, which is needed to specify the size of T. 
However, if (i[k],k G A' corresponding to the trigonometric 
polynomial p is a solution to whose support set A' is 
strictly smaller than the assumed support set A, then any 
multiple u = p • y having coefficients e[k] = {d ^ g)[k\ 
supported within A, is also a solution. This implies that if 
we consider a larger filter size than required by the minimal 

















































polynomial, i.e. more columns in T than the number of 
coefficients in d, the matrix T will be low-rank. 

The preceding discussion suggests we may pose the recov¬ 
ery of the signal as a structured low-rank matrix completion 
problem, entirely in the Fourier domain: 

/ = argmin rank(T[^]) subject to ^[k] = /[k],Vk G F 

9 

( 18 ) 

We note this approach is still “off-the-grid” in the sense that 
we may recover a discrete image at any desired resolution by 
extrapolating / to this resolution in Fourier space, and apply¬ 
ing an inverse DFT. To address the case of noisy measurements 
and model mismatch we propose solving the convex relaxation 
of ([T^: 

/ = argimn ||T[5]||* + A||Pr(5 -/)||2 (19) 

9 

where || • ||* denotes the nuclear norm, i.e. the absolute sum 
of singular values, A is a tunable parameter, and Pr is the 
projection onto the sampling set F. A standard approach to 
solving © is by an iterative singular value soft-thresholding 
algorithm, which requires an SVD of the estimate T\g] at 
each step. Due to the size of T\g], such an algorithm is 
computationally prohibitive in this case. Instead we use the 
SVD-free algorithm proposed in |T^ , which involves intro¬ 
ducing auxiliary variables U G and V G via 

the well-known relation ||X||* = minx=uv^ l|U|||^ + II^IIf^ 
and enforcing the constraint T[^] = UV^, with the ADMM 
algorithm. 

In Fig. we demonstrate the ability of the proposed 
algorithm to recover a piecewise constant signal from few of 
its uniform low-resolution Fourier samples. We experiment 
on simulated data obtained from analytical MRI phantoms 
derived in |T3| . We extrapolate from 65 x 49 = 3185 analytical 
Fourier samples of the Shepp-Logan phantom to a 256 x 256 
grid (^20-fold undersampling), and recover the signal by 
performing a inverse DFT. Note that the ringing artifacts 
observed in the recovery are to be expected due to fact we are 
recovering exact Fourier coefficients of the signal, and could 
be removed with mild post-processing. 

VI. Conclusion 

We propose an extension of the annihilating filter method to 
a wide class of 2-D piecewise smooth functions whose edges 
are supported on level set of a band-limited function. This 
enables us to recover an exact continous domain representation 
of the edge set from few low-frequency Fourier samples. In 
the case of piecewise constant signals, we derive conditions 
of the necessary and sufficient number of Fourier samples 
to ensure exact recovery of the edge set. Lastly, we prosed 
one-stage algorithm to recover piecewise constant signals by 
extrapolating the signal in Fourier domain. We demonstrate 
that we may accurately recover MRI phantoms from few of 
their low-resolution Fourier samples. 



(a) Fully sampled (b) Zero-padded (c) Low-rank 



Fig. 2. Recovery of Shepp-Logan phantom on a 256x256 grid from 
65x49=3185 Fourier samples (?^20 fold undersampling). The top row 
shows the spatial domain images, while the bottom row shows the 
Fourier transforms of the images (log scale). 
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